%----------------
% makeruns : generate two fake data sets for testing the routine above.

function [a, b] = makeruns(testmode, pflag)

start = 100;
baseper = 20; % 50 Hz.
stimdur = 1000;
ntestrep = 50;
a=[];
b=[];
switch(testmode)
    case 1 % exactly synchronized, and regular
        for i = 1:ntestrep
            a{i}.latency = [start:baseper:stimdur];
            b{i}.latency = [start:baseper:stimdur];
        end;

    case 2 % one is regular, the other has random jitter
        for i = 1:ntestrep
            X = [start:baseper:stimdur];
            a{i}.latency = X;
            b{i}.latency = X + 5*randn(size(X)); % normal random distrubiton on b dataset
        end;

    case 3 % both are regular, but both have independent random jitter
        for i = 1:ntestrep
            X = [start:baseper:stimdur];
            a{i}.latency = X + 5*randn(size(X));
            b{i}.latency = X + 5*randn(size(X)); % normal random distrubiton on b dataset
        end;

    case 4 % both are random and independent
        for i = 1:ntestrep
            X=1+randpois(baseper*3, floor((stimdur-start)/(baseper*3))); % 1 adds fake refractory period
            a{i}.latency = start + cumsum(X); % corresponds to about 50 s/s
            X=1+randpois(baseper*3, floor((stimdur-start)/(baseper*3))); % separate process
            b{i}.latency = start + cumsum(X); % corresponds to about 50 s/s

        end;
        ddur = 100;

    otherwise
        return;
end;

if(pflag) % then try plotting the spike trains.
    h = findobj('tag', 'makeruns');
    if(isempty(h))
        h=figure('tag', 'makeruns');
    else
        figure(h);
        clf;
    end;
    for i = 1:ntestrep
        ya=i*ones(length(a{i}.latency));
        plot(a{i}.latency, ya, 'rd');
        hold on;
         yb=i*ones(length(b{i}.latency));
        plot(b{i}.latency, yb, 'b+');
    end;
    set(gca, 'Xlim', [0 1000]);
    set(gca, 'Ylim', [0 ntestrep+1]);
end;



% RANDPOIS(p,n)
%   This function creates a random variable with a Poisson distribution
%   with parameter 'a'.  If a second argument is used, a vector of
%   'n' poisson random variables is created.
%
% See also RAND, RANDN, RANDUNIFC, RANDEXPO, RANDGEO, RANDGAUSS

function out = randpois(a,n)

if nargin == 1
    i = 0;
    f = exp(-a);
    p = exp(-a);
    u = rand;
    while f<=u
        p = p * (a / (i + 1));
        f = f + p;
        i = i + 1;
    end
    out = i;
end

if nargin == 2
    for k=1:n
        i = 0;
        f = exp(-a);
        p = exp(-a);
        u = rand;
        while f<=u
            p = p * (a / (i + 1));
            f = f + p;
            i = i + 1;
        end
        out(k) = i;
    end
end